Numerical and experimental investigation of multi-species bacterial co-aggregation

This paper deals with the mathematical modeling of bacterial co-aggregation and its numerical implementation in a FEM framework. Since the concept of co-aggregation refers to the physical binding between cells of different microbial species, a system composed of two species is considered in the modeling framework. The extension of the model to an arbitrary number of species is straightforward. In addition to two-species (multi-species growth) dynamics, the transport of a nutritional substance and the extent of co-aggregation are introduced into the model as the third and fourth primary variables. A phase-field modeling approach is employed to describe the co-aggregation between the two species. The mathematical model is three-dimensional and fully based on the continuum description of the problem without any need for discrete agents which are the key elements of the individual-based modeling approach. It is shown that the use of a phase-field-based model is equivalent to a particular form of classical diffusion-reaction systems. Unlike the so-called mixture models, the evolution of each component of the multi-species system is captured thanks to the inherent capability of phase-field modeling in treating systems consisting of distinct multi-phases. The details of numerical implementation in a FEM framework are also presented. Indeed, a new multi-field user element is developed and implemented in ANSYS for this multiphysics problem. Predictions of the model are compared with the experimental observations. By that, the versatility and applicability of the model and the numerical tool are well established.

www.nature.com/scientificreports/ From a clinical perspective, the development and uncontrolled expansion of oral biofilms can lead to diverse pathologies, primarily local but potentially with subsequent systemic complications. Oral biofilms are responsible for the most prevalent infectious disease worldwide, dental caries and periodontal disease, and the most common implant-associated complications, peri-implant diseases. Odontogenic complications include pathologies at diverse anatomic locations including the cardiovascular system, gut, respiratory tract, reproductive tract, or even brain. Physical contact between biofilm members contributes to biofilm structure that facilitates higher resistance to immune response and antimicrobial treatment. On the one hand, co-aggregation can directly contribute to virulence in a process of synergistic invasion 12 or invader recruitment 13 . On the second hand, receptors involved in co-aggregation can be used to target pathogens, with the use of antimicrobials and probiotics 14 . Given the critical role of cell-cell contact in biofilm development and maturation, it can be speculated that a better understanding of this process could result in improved prophylaxis and anti-biofilm therapy. Yet, the holistic approach to oral microbial ecology is challenged by the high diversity and complexity of microbial communities. High-resolution numerical simulations can help to overcome current methodological limitations and improve a system-level understanding of oral biofilms.
The mathematical modeling and numerical simulation of bacterial colonies (biofilms) have drawn a great deal of attention among researchers over the last two decades. In particular, the spatiotemporal evolution of bacterial aggregates is of great interest in the computational biomechanics community. A common and favorable modeling scheme is Individual-based modeling (IBM) which is very flexible in terms of accommodating complex interactions between the discrete agents as well as sophisticated biological, chemical, and physical processes associated with a single cell such as microbial cell division, decay, motility, etc. Moreover, this approach is readily scalable in the sense that an extremely large number of microbial populations can be handled using parallel solvers and domain decomposition due to the explicit nature of their time integration scheme. As the outcome of such endeavors, one can name several open-source computer codes such as iDynoMiCS 15 and Simbiotics 16 , NUFEB 17 and also MIRODIM 18,19 which are capable of simulating the development of microorganism in a multidimensional and multi-species fashion. IBMs are normally of hybrid type meaning that besides the bacterial species, which are treated as discrete agents, the nutrient availability is described using a continuum approach. Along with IBMs, fully continuum-based modeling of biofilms is also mathematical modeling and the numerical simulation of the multi-species biofilms have not yet been sufficiently considered in the literature. Although most species that constitute dental plaque are known individually, their collective behavior and group interactions are still under investigation. In order to understand the spatiotemporal structure of such a bacterial community one has to identify and understand the interaction, partnership and communication among the members of this community. In particular, the mathematical modeling of such phenomena is really sparse in the literature.
When it comes to mathematical modeling, the growth of biofilms can be viewed from two different perspectives. In the first approach, the growth is described as volume expansion and consequently, a "mechanical movement" of the boundary of the growing matter. This "explicit" movement of the biofilm can be captured, in the simplest fashion, using a potential flow which is isotropic, see for example 20 and 21 . Alternatively, in a more complex framework, one can use anisotropic kinematics that governs the deformation of the biofilm as a deformable viscoelastic solid instead of a fluid, see 22 and 23 . This notion of growth has also been employed to model various growth phenomena in soft tissues other than bacterial aggregates. Examples are tumor growth 24 , and aneurysm as a result of abnormal arterial growth 25 . In the second approach, the growth of biofilm can be characterized using the "propagation" of a phase-field variable. A classical example is the classical diffusionreaction equation in which the propagation of the wavefront can be interpreted as the growth of the colony of species, see 26,27,[28][29][30] . Another example is Fisher's equation 31 which is, indeed, a transient diffusion-reaction model for predicting population dynamics and pattern formation. One can generalize this framework by using a set of diffusion-reaction equations in order to account for different reciprocal interactions, such as symbiosis or competition, between multiple species. For instance, the diffusional Lotka-Volterra system 32 is associated with a set of PDEs capable of describing the interaction between two species which are, in fact, prey and predator.
In this work, we use the second approach in order to mathematically describe a system composed of two different bacterial species which not only grow but also co-aggregate in the presence of each other. Moreover, it is shown that this approach is exactly equivalent to employing a so-called "phase-field" model for capturing the behavior of such biological systems. The mathematical model is implemented using a fully implicit and monolithic scheme within a FEM framework. A multi-physics user element is developed from scratch to be invoked by the ANSYS solver. Several numerical examples are presented to show the robustness of the proposed numerical tool in handling multi-species biological systems.

Mathematical modeling of the co-aggregation of two bacterial species
From a biological point of view, co-aggregation is a crucial process in the formation of multi-component biofilms. Genetically distinct bacteria adhere to each other, ' co-aggregate' , and can replicate faster, see 33 . Due to co-aggregation, the metabolism of the species is altered compared to the monocultures, see 3 . In order to mathematically describe the development of biofilm under the influence of co-aggregation, the basic assumptions are: • Two species of bacteria are taken into consideration. One can readily extend the model to multiple species without any restriction. • For the sake of simplicity, only a single nutrient is assumed to nourish both species. The nutrient is, indeed, dissolved in the fluid surrounding the biofilm. One can extend the model to multiple nutrients. • Since the fluid flow is not modeled explicitly, the transport of the nutrient is governed by the diffusion mechanism. This is in line with a very small flow velocitythat is either induced by the biofilm growth (in the order of a few microns per day) or governed by quasi-hydrostatic conditions. As a result, the advection mechanism www.nature.com/scientificreports/ for the transport is neglected. Such a situation can occur in specific anatomic locations, e.g., deep pockets around implants or teeth, where salivary flow is strongly reduced.a very small Peclet number (Pe) allows such simplification regarding the transport mechanism. One can estimate that Pe = Advection rate Diffusion rate = Lu D → 0 in which L, u and D refer to the characteristic length, characteristic fluid velocity, and nutrient diffusivity, respectively.
• The generation of so-called byproducts, as a result of the metabolic activities of the bacteria, is not modeled explicitly. It is known that metabolic products, such as lactic acid, can have either a negative or positive influence on the growth of different components in multi-species systems 34 . • The model of the co-aggregation is mathematically based on the "proximity" of species. Irrespective of the underlying biological mechanism, the model assumes that bacterial growth is boosted if particular species are "sufficiently close" to each other. In this sense, co-aggregation can be perceived as cooperation. • It is postulated that the growth of a biofilm colony is in the direction of the maximum nutrient gradient. It will be clarified more in the mathematical description of the proposed phase-field model. • The initial seeds of species are randomly generated in the geometrical model. The model does not capture the very early colonization of the surface by the planktonic bacteria. Such a process needs to be modeled using agent-based methods with regard to the length scale of bacteria and the surface roughness, see 35 .
The mathematical model is based on these assumptions. For the sake of simplicity, the schematic view of the computational domain is considered to be two-dimensional. The variables φ 1 and φ 2 are taken to represent species 1 and 2 respectively. These two variables can be perceived as the density of bacteria. If they are equal to zero at a certain location, it implies the absence of that species. Finally, the sharp interface between zero-valued and one-valued regions characterizes the boundary of the species. Figure 1 depicts the regions where φ 1 and φ 2 change spatially.
One can see that both species may coexist (both φ 1 and φ 2 amount to 1) in certain regions. Such regions are exactly the aggregation zones. The aggregation is captured using a new variable called α . This variable can be evaluated by applying a Boolean operation of "AND" type on φ 1 and φ 2 . This variable is treated like a phase-field variable. It means that α = 1 corresponds to the regions in which both φ 1 and φ 2 are present and consequently the co-aggregation can take place.
Finally, the concentration of the nutrient is denoted by the variable c that obeys a diffusion-reaction equation. It is obvious that the regions occupied by one of the species, namely φ 1 and/or φ 2 equal to 1, play the role of sink term (reaction) in the diffusion-reaction equation. Furthermore, the values of nutrients on the boundaries can be prescribed/controlled in the context of boundary conditions. Phase-field type equations for species and co-aggregation. As mentioned before, the variable α is captured using a phase-field equation. The well-known Allen-Cahn 36 type of phase-field equation can be written as in which the generic variable can be replaced by either φ 1 , φ 2 or α . The dot over in the first term signifies the time derivative. The parameter ε regularizes the sharp change of φ between the two phases. f ′ ( ) is a short notation for ∂f ( ) ∂ . The function f ( ) is a so-called double-well potential and characterizes the barrier that has to be overcome in order to have a phase transformation (here, phase transformation means a region that used to be free of species, namely = 0 , is colonized by that species, namely = 1 ). The description of each term in the equation is self-explanatory and standard in any application of phase-field modeling.
A common choice for the double well potential is Driver (source) of phase-field , Figure 1. Phase-field approach for a multi-species system, φ 1 and φ 2 represent species 1 and species 2 while α refers to the congregation state. www.nature.com/scientificreports/ in which M is proportional to the local maximum of the function at φ = 1 2 between the two wells (minima) at = 1 and = 0 , see Fig. 2. By degenerating the problem to 1D one can better understand the idea behind the phase-field approach. Neglecting the time-dependent as well as the source term, in 1D, one can realize that the so- (1). Here l = ε 2 M characterizes clearly the "interface width". It captures a smooth but sharp transition between the two phases. In a similar fashion, one can show that the "interface energy" is proportional to √ ε 2 M . In other words, it represents the resistance against phase change. If the driving force (the last term in Eq. 1) overcomes this threshold, the phase change takes place. Actually, ε and M are characteristic parameters that govern the underlying physical process pertaining to the phase-field equation. The readers are referred to 37 for more elaboration on the structure of the phase-field equation.
As stated before, the co-aggregation ( α ) can be captured using the phase-field equation. It is noteworthy to mention that this type of phase-field modeling is, in fact, a transient diffusion-reaction equation. To restore the classical diffusion-reaction equation one can eliminate the term concerning double well potential (by setting M = 0 ). In fact, the evolution of the two species' density is governed by the classical diffusion-reaction equation.
To do so, one can substitute for φ 1 or φ 2 and set M = 0.
The crucial part of the phase-field equation is the definition of the source terms for each field. The source terms corresponding to the four variables φ 1 , φ 2 , α, c are denoted by , respectively. That gives rise to the evolution of the system. Its mathematical representation needs to be physically meaningful. As stated in the assumptions, the "physically meaningful rationale" is that the expansion of bacterial colonies is preferred in the direction of maximum nutrient change, namely nutrient gradient. Nutrient metabolism provides energy and building blocks for all cell replications. Consequently, in locations where access to nutrients is hindered, replication is slower or stops completely. In the following, appropriate source terms are considered.
The function S 1 (φ 1 , α, c) can be proposed as where R s1 is a parameter that controls the magnitude of the source term and τ characterizes the degree of growth promotion due to coaggregation. It is hypothesized that the growth (bacterial replication) is boosted if α is evolved (co-aggregation takes place). As long as different species have positive effects on each others' fitness/ growth, closer proximity introduced by co-aggregation can facilitate stronger synergy. The function H(φ 1 − φ cri ) is the Heaviside step function with a jump at critical value φ cri . A general definition of H(· − · cri ) with a jump at · cri is One can interpret the parameter φ cri as the minimum level of bacterial density that prevents the bacterial colonies from vanishing, in the long run. The source term, which can be interpreted as bacterial replication, is proportional to bacterial density and nutrient availability. The source term for species 2 is similar to species 1 The source term for the co-aggregation variable α can be written as  www.nature.com/scientificreports/ in which the term ∇c |∇c| is a unit vector pointing to the direction of maximum change in the nutrient. This is, indeed, the preferred direction for the expansion of colonies. Taking a closer look at the structure of S α (φ 1 , φ 2 , α, c) , one can notice that it serves as an advection term in the form of v · ∇α 1 and v = R α ∇c |∇c| being the advection velocity.
It means that the boundary of the bacterial colonies at the co-aggregation zone is advected in the direction of v.
For the nutrient c, the source term should be negative, because the nutrient is consumed.
in which the parameter R c is the consumption parameter. Obviously, the transport of the nutrient is governed by the standard and classical time-independent diffusion-reaction equation ( M = 0 ). For the nutrient, the parameter ε 2 in Eq. (1) is, in fact, the nutrient diffusivity. To make a distinction, we replace ε 2 with D n and D b in the case of nutrient and bacterial species, respectively.

Method
All methods were performed in accordance with the relevant guidelines and regulations.
Numerical implementation using FEM. Since the adopted numerical method is based on the standard Galerkin FEM, one needs to convert the differential equations into their respective weak forms. One can express the weak form of four governing equations for the existing four unknowns ( φ 1 , φ 2 , α and c) as follows: One should note that variation of the variables ( δφ 1 , δφ 2 , δα, δc ) serve as the test function in the context of variational approach. Furthermore, (·) acting between the variable gradients is indeed the inner product. In line with standard FEM, the objective is to find the field variables with so-called C 0 continuity . Upon applying integration by parts to the integrals in Eqs. (8)- (11), the boundary terms defined on the boundary ∂B emerge. Furthermore, a backward (implicit) Euler scheme is used for the time discretization which replaces the time derivatives. This yields where N is the normal vector of the boundary surface. It is an acceptable assumption that the boundary flux term for the phase-field is assumed to be zero. Moreover, D∇c · N refers to the nutrient flux at the boundary. The superscript ( n − 1 ), which appears due to time discretization, refers to the values of variables at the previous time step. The implementation of the multi-field (species 1, species 2, co-aggregation, and nutrient concentration) threedimensional problem at hand was carried out using AceGen, see 38 , which is based on automatic differentiation (hybrid symbolic/numeric differentiation). Once the weak forms are linearized, the stiffness matrix is extracted and tailored to a user element in the FORTRAN language that can be implemented in any FEM code using standard discretization schemes. The numerical test cases are performed on two-dimensional models for the sake of computational efficiency. Due to the nature of phase-field equations, one needs to use sufficiently fine mesh. Even in two-dimensional examples, a mesh size of 0.5 µm leads to approximately 2 × 10 4 nodes which means around

B
[D n ∇ · (∇c)δc + S c (φ 1 , φ 2 , c)δc +ċ δc] dv = 0. www.nature.com/scientificreports/ 8 × 10 4 degrees of freedom (DoF) in the final model. One can easily understand that for a three-dimensional model with the same mesh size, the total number of DoFs exceeds 10 7 and is computationally prohibitive. Each node of the FEM mesh has four degrees of freedom all of which are scalar-valued. The first and second ones are dedicated to φ 1 and φ 2 representing the density of species 1 and species 2, respectively. The third is α which captures the co-aggregation degree of the two species. The fourth one is allocated to the concentration field c. The solution of the global system yields the updated value of the field variables, with the assumption that the previous values ( φ n−1 1 ,φ n−1 2 ,α n−1 and c n−1 ) are known. In the next section, several numerical examples are presented to see the applicability of the numerical tool. The post-processing and visualization of the results were conducted in Paraview.

Numerical and experimental test cases.
For the numerical simulations, all geometrical data as well as material parameters are summarized in Table 1. the temporal unit (denoted by "Time" in the table) is nondimensionalized using growth characteristic time scale which is in the order of a few hours, see Fig. 3. It means that, from a numerical point of view, "Time" varies between 0 and 1. Test case 1: single colony, single species. In the first example, different patterns of growth are studied.

Approval for human/animal experiments.
A single tiny colony is assumed to be located at the center of the domain. The initial nutrient distribution is uniform and constant. The numerical investigation shows that the ratio between the nutrient and biofilm diffusivity ( r = D n D b ) plays an important role in the morphology of the grown biofilm. Hereafter, two distinct cases are investigated. It is assumed that the entire domain is enriched by a uniform initial nutrient. There is no external provision of the nutrient. Hence, the nutrient can only be consumed. To validate the reliability of the mathematical model, experimental observations are also conducted.
Simulation: diffusion-limited growth (dendritic patterns). If r D is relatively small, the biofilm tends to branch and form a dendritic shape. This is referred to as diffusion-limited growth in the literature. The formation of a finger-like structure, from a mathematical point of view, is rooted in the fluctuation-induced interface instability, see 39 . Such patterns are also observed in other phenomena such as multiphase flows 40 and solidification process 41 . The formation of all irregular morphology with lobate margins can be explained in the context of nutrient absorption. It means that as a result of the scarcity of nutrient flux due to low diffusion through the biofilm, the bacteria increase the absorption surface via forming lobes. One can see that the nutrient is predominantly consumed at the locations where the biofilm exists and the far-field is not influenced by the growth. This is due to relatively low nutrient diffusivity that can not compensate for the consumed material. Looking at the www.nature.com/scientificreports/ pattern of nutrient consumption, one can realize that the nutrient is only depleted at the regions where biofilm exists. Figure 4 shows the numerical simulation of this scenario that is in qualitatively good agreement with the experimental observations reflected in Fig. 3

part b.
Simulation: diffusion-unlimited growth (compact pattern). If the diffusion of the nutrient is fast, the biofilm grows in a compact and spherical fashion, see Fig. 5. Unlike the previous example, the fast diffusion of nutrients ensures optimal supplies for the bacteria. Hence, the overall shape of the bacterial colony remains circular with entire margins. Here, due to the high diffusivity, the nutrient is transported from the whole domain to the biofilm-colonized regions. This leads to the gradual depletion of the nutrient in the whole domain, unlike the previous case. Experimental investigations also comply qualitatively with the mathematical model, see Fig. 3 part a.

Test case 2: multiple colonies, two species, with interaction (co-aggregation). Experiments:
sample preparation for co-aggregation and biofilm formation. The bacteria were transferred to a co-aggregation buffer (1 mM Tris, 150 mM NaCl, 0.1 mM CaCl2, 0.1 mM MgCl2) and the optical density at 600 nm was set to www.nature.com/scientificreports/ 1.0. That is equal to a cell density of about 109 cells per ml. The co-culture inoculum was generated by mixing the same volumes of two species and resulted in the formation of co-aggregates. Monoculture and co-culture inocula were vortexed for 10 seconds, diluted ten times in co-aggregation buffer followed by mixing again. Then 3 ml of suspension was placed in 6-well plates (Thermo Fischer Greiner Bio-one Polystyrene 6-well Multiplates) and allowed to coat for 1 hour. Following this, the non-attached bacteria were removed from the wells and 3 ml of THBY medium was added. The bacteria were allowed to grow for 8 hours from this time point. After 8 hours of growth biofilms were stained with the LIVE/DEAD BacLight Bacterial Viability Kit (Invitrogen); rinsed, fixed in 2.5% glutaraldehyde solution in phosphate-buffered saline, and examined using CLSM (Leica TCS SP2, Leica Microsystems, Mannheim, Germany). For the imaging purpose, 488 nm laser was used with an emission range of 500 -545 nm for SYTO9 and 590 -680 nm for propidium iodide (PI). For the biofilm in each well, five representative images were acquired with optical sections of 3 µm . The Imaris Cell Imaging Software package (Imaris x64, 6.2.1, Bitplane AG, Zürich, Switzerland) was used to analyze images. A three-dimensional reconstruction smart view was generated, with green (SYTO9), red (PI) and yellow (co-localized SYTO9 and PI) fractions. Non-permeable cells were stained in green. Yellow and red signals were considered to represent permeable cells www.nature.com/scientificreports/ since PI penetrated the membrane of these cells. The experiment was performed twice. Each time for each group a total of two wells and ten images were generated. The experimental observations are illustrated in Fig. 6.
Simulation. In this test case, unlike test case 1, multiple colonies of two species are considered. The co-aggregation mechanism also exists and contributes to growth boosting. In order to replicate numerically the experimental observations shown in Fig. 6, the entire computational domain is randomly filled with initial colonies of both species in a uniform way. Then, we let the system evolve in time. The formation of the bacterial cluster in the regions with high proximity of the two species is observed as a result of co-aggregation. Figures 7 and 8 show snapshots of both species being represented by tiny colonies with different colors (gray and cyan represent species 1 and species 2, respectively). Additionally, the regions where co-aggregation has taken place are pictured using the red rods. These two figures show the development of the colonies over the course of time for two different conditions discussed in the previous example, namely diffusion-limited and -unlimited www.nature.com/scientificreports/ growth. Similar to the previous test case, the appearance of the dendritic structure is the result of a diffusionlimited scenario (small r D ), while round and compact islands of colonies are a sign of diffusion-unlimited growth. These figures also demonstrate the nutrient distribution in the course of time. One can see that the nutrient consumption map is significantly different in the two cases. All in all, the numerical predictions show plausible qualitative compliance with the experimental observations. The average nutrient consumption is reflected in Fig. 9. While the nutrient concentration drops gradually in the diffusion-limited case, it is rapidly depleted in the diffusion-unlimited condition. In both cases, the profile of the average nutrient is monotonically decreasing because the uniform initial concentration is only consumed by the biofilm. It is also reflected in the average density of the biofilm that monotonically increases in both cases, however with different rates.

Test case 3: growth under the directional provision of nutrient with/without reciprocal interaction (co-aggregation).
This example is identical to the previous one with the difference that the nutrient is supplied constantly from a particular location (the top side, here). It leads to a prescribed gradient profile of the nutrient and consequently a preferred direction of colony expansion. To quantify the growth process, one can look at the average density of the biomass produced in time. It reflects how the available empty space is occupied by the two-species biofilm.
The experimental setup for such a test case is interesting and quite elaborate because the continuous and purely diffusion-driven feed of an aqueous environment from a specific location (one side) needs special care and measures. Any perturbation in the solution fluid can result in convective movement and mixing phenomenon which destroys the intended gradient in the nutrient availability. Unlike the experimental setting, the numerical modeling of the problem is free of such practical complications. One needs to apply the desired boundary condition on the intended regions. Hence, the experimental investigation is ruled out for this test case.
The numerical simulation shows how the nutrient gradient affects the growth pattern of the biofilm. In fact, the final morphology is a result of both co-aggregation and nutrient availability. The biological activity of the bacteria is intensified according to the availability of co-aggregation partners as well as nutrient. Figs. 11 and 12 demonstrate snapshots of bacterial development in a dual-species system. The regions in which both species are coexistent are the dominant sites of bacterial growth which contribute to transforming tiny colonies into larger  www.nature.com/scientificreports/ structured clusters of bacteria (biofilm). The clusters tend to form in the regions with high amounts of nutrients, namely the regions in the proximity of feeding sites at the top. Apart from nutrient feed, comparing Figs. 11 and 12 reveals how substantially the co-aggregation contributes to the rapid biofilm development. In Fig. 10, the profile of the average nutrient and biofilm density is plotted with respect to the time. Due to the continuous provision of the nutrient and tiny initial colonies, the average value of the nutrient concentration increases first (primary ascending phase) and then decreases to a steady state constant value (secondary descending phase). The steady-state condition is reached when the rate of nutrient feeding is balanced with the rate of nutrient consumption when sufficient biofilm is developed. In the absence of co-aggregation, the development of biofilm, which determines the consumption capacity, is significantly slower. Consequently, the  www.nature.com/scientificreports/ nutrient feed outweighs the nutrient consumption leading to an increase in the average nutrient and a delayed succeeding decrease.

Conclusion
A novel mathematical model based on the phase-field method is proposed for the simulation of multi-species biofilms. The co-aggregation of the species is taken into account in a system of multi-species colonies. The complex morphology of multi-species colonies observed in the experiments can be reproduced using the presented mathematical model. The formation of such complex structures is driven, on one hand, by "nutrient availability" and, on the other hand, by "species proximity". Since the proposed model is based on the phase-field approach, one does not need to track the moving interfaces of the multi-species and multi-colony biofilms. Handling multiple and moving boundaries is the main superiority of phase-field modeling in this application. Additionally, www.nature.com/scientificreports/ the merging of multiple colonies is treated naturally. It turns out that phase-field modeling is a successful and applicable approach for this biological problem at hand. The nature of biofilm and the underlying processes are so complex that probably four continuous isotropic scalar variables ( φ 1 , φ 2 , c and α ) are not sufficient to capture the real bacterial morphology to a good extent. This is why quantitative validation of the numerical results against the experimental observations in this work is a challenging task. Nevertheless, mathematical modeling, despite simplifications, can help us to understand the "fundamental underlying mechanisms" that are involved.
This work can be extended and improved in several directions. A system composed of more than two species is of great interest and can be modeled just by adding phase-field equations. Furthermore, the interaction between the species can be generalized to not only reciprocal dual interaction but also to multilateral metabolic communications. Incorporating anisotropic growth function can also enhance the model. Finally, an efficient parallelized numerical architecture (such as MPI) can be employed to solve the model in 3D. These extensions are left for future work.